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Abstract 

The non-relativistic Schrodinger equation with the linear and Coulomb potentials is solved 
numerically in configuration space using the relaxation method. The numerical method 
presented in this paper is a plain explicit Schrodinger solver which is conceptually simple 
and is suitable for advanced undergraduate research. 



1 



1 Introduction 



This work came out of an undergraduate research project at the University of Wisconsin- 
Milwaukee. The analytical solution of the hydrogen atom is a typical topic in advanced 
undergraduate or beginning graduate quantum mechanics which utilizes the series expan- 
sion technique. The analytic solution does not only illustrate the mathematical approach 
of solving differential equations, it also provides a benchmark for testing numerical meth- 
ods. Initially the numerical solution of the hydrogen atom was intended to offer the 
undergraduate students in our department (the second and third authors of this paper) 
an opportunity to learn the elements of scientific computation and basic research strate- 
gies. The students began by learning the theory of partial differential equations, linux 
programming in C, and basic numerical methods. In order to give the students a taste 
of original theoretical research, they were asked to check the eigenvalues generated by 
a Nystrom momentum space code from a new paper with the configuration space code 
which they helped to develop. This paper documents the r-space code for the numerical 
solution of the non-relativistic Schrodinger equation (NRSE) with the Coulomb and linear 
potentials. It is hoped that this work is helpful to those students and teachers who want 
to integrate numerical methods with a traditional quantum mechanics curriculum. 



2 Schrodinger Equation in configuration Space 



The basis of the wavefunction of the Schodinger equation in configuration space is taken 
to be 

(j)imir) = ^YiU^) (1) 



The radial part of the equation for a hydrogen atom is |TI| 
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A computer cannot integrate r from zero to infinity. Therefore we must map [0, oo) — >• 
[0, 1] by 

r = (3) 



It implies that 
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From now on, a new symbol for the radial wavefunction will be used, i.e. y = Ri. 
To transform the Schrodinger equation to the new space, we first transform the second 
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derivative, 
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By substituting Eq. H] into Eq. H] and letting h ^ 1 and 
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/i = 0.5107208 X 10*^ el/, 



in natural units, the Sclirodinger equation is transformed as 
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y = 0. (7) 



With the benefit of hindsight, we know that the non-zero portion of the solution of Eq. 0] 
will be tightly bound to a narrow margin near x ~ 0. In order to improve the numerical 
accuracy and the quality of the plots, we like to focus on the non-zero portion of the 
solutions. A new configuration variable z and a redefinition of x 
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are used to zoom into the small x region. Eq. [0 is modified as 



d^y 2 dy 1 



dx"^ 1 — X dx (1 — xY 



2nal[E + 



1 — X e 



X ao 



1 — X 



X 



1(1 + 1] 



(9) 



y = 0. (10) 



3 Relaxation Method 

The shooting method typically shoots from one boundary point to another using Runga- 
Kutta integration. In the case of the Schrodinger equation, the boundary points at x = 
and X = 1 are both singular. A one-point shoot will not converge when the code marches 
toward a singularity. A two-point shoot marches from both boundary points at a; = 
and X = 1 to match a third boundary point somwhere in between. Even if it works, 
two-point shoot requires too much a priori knowledge of the wavefunction. Furthermore 
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any Runga-Kutta based alogrithms, such as the shooting methods, will fail under normal 
circumstances because either (1) the integration tends to blow up to infinity when shooting 
from x = or (2) the integration is identically zero when shooting from x = 1, given 
the boundary conditions y{l) = y'{l) = 0. An embedded exponentially- fitted Runga- 
Kutta method is reported to work. But it is beyond the scope of an undergraduate 
research project. The relaxation method, on the other hand, will in principle handle two 
singular boundary points. In this paper, we adapt the relaxation codes used in Numerical 
Recipes In order to establish a continuity between this paper and the said reference, 
the same notations are used. 

The basic idea of the relaxation method is to convert the differential equation into 
a finite difference equation. Error functions are then considered. For instance, the er- 
ror function at the A;-th mesh point of the j-th differential equation (one of coupled 
differential equations) 

^ = 9ixk,yu...,yN) (11) 

is simply 

Ej,k = yj,k - yj,k-i - {xk - Xk-i)g{xk,Xk-i,yi,k,yi,k-i, ■ ■ ■ ,yN,k,yN,y-i)- (12) 

The difference of Ej ^. is approximated by a first order Taylor series expansion, such that 

N 2N 

Si,n^yn,k-l + X! Si,n^yn~N,k = ~Ei^k, (13) 
n=l n=N+l 

where 

'Ji,n — T{ ) iJi,n+N — ■ [i-^) 

oyn,k~i oyn,k 

The crux of this paper is to show how the matrix elements Si^k of the non-relativistic 
Schrodinger equation are calculated and how to streamline the relaxation codes. The 
ranges of the sums in Eq [^] are split over [1, N] and [A^+1, 2N] because of the peculiarity 



of the relaxation codes used in Numerical Recipes. We will not delve into the details of 
the codes in this paper. Readers are encouraged to refer to Numerical Recipes for a full 
explanation. As a motivation, it suffices to say that the main idea of the relaxation method 
is to begin with initial guesses of yj^k and relax them to the approximately true values 
by calculating the errors Ei^k to correct yj^k iteratively. yj^k are components of a solution 
vector in an AM-dimensional vector space. Intuitively we can think of the relaxation 
process as rotating an initial vector into another vector under the constraints of Ei^^- 
Since y = is a trivial solution, the relaxation process has a tendency to diminish those 
components of the solution vector that correspond to the wavefunction and its derivative. 
We will return to this point later in the paper. 
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Let M — 1 be the number of mesh points and h is step size, then 

1 



h 
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We define the vector yi as 



{k-l)h, 
yixk). 



yi = y 
y2 = y' 
. y3 = E 



(15) 



(16) 
(17) 



such that 



Vi 



y'2 



y2 



2 1 

y2 - 



1 — X 



{i-xY 



2^ h/3 + 



1 — X e 



X ao 



1 — X 



X 



yi 







A first order finite derivative at Xk is given as 



/ yk - yk-i 



In the relaxation formahsm, Ei^^ are to be zeroed. In this case, it is easy to see 
Eq. can be rewritten as 
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The corresponding Sij S'-matrix elements are 
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Similarly, Eq. is rewritten as 
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At last, Eq. |^ is rewritten as 
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and the associating S^j matrix elements are 
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At the first boundary point x = 0, we have rii = 1. The boundary condition is y{0) = 0. 
It translates to 

£^3,1 = 1/1,1. (43) 

The only non-zero S'-matrix element is 

^3,4 = 1. (44) 

At the second boundary point a; = 1, we have 77-2 = 2. The boundary conditions are 
y(l) = y'{l) = 0, or equivalently 



Ei^M 
E2,M 



The only non-zero matrix elements are 



>S'l,4 
'S'2,5 



yi,M, 
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1. 
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The hydrogen atom finite difference equation can be generalized to any Coulomb 
pentential by substituting — > Ac- To adapt the code to solve NRSE with a hnear 
potential, the only modifications needed are 



in Eq. p 



1 - 



X 
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(49) 



(50) 



in Eqs. [29,30|, and the omission of all uq from the matrix elements. 
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4 Exact Solution of the Hydrogen Atom 

The non-relativistic hydrogen atom is solved exactly by series expansion. The first few 
radial wavefunctions are found to be 

^io(r) = ^7^e"^^^ (51) 

R2o{r) = -^4r{l-lir)e-^^\ (52) 

i?2i(r) = -^4r^e-^^\ (53) 

where 

ln = ^ (54) 

With the redefinition of variable by Eqs. HO] and the omission of normalization constants, 



Eqs. 51 -53 1 can be rewritten as 



i?io(^) ~ ^e-^ (55) 
R2o{z) ~ 2(l-2)e-t, (56) 
R2i{z) ~ z'e-i. (57) 



Eqs. []55|-p7| will be used to test the relaxation codes in Section 0. 

5 Exact /S-state Solution for the Linear Potential 

The S'-state eigenvalue of the non-relativistic Schrodinger equation with a linear potential 
can be solved exactly in configuation space. We shall use the analytic results to check our 
numerical results. The non-relativistic Schrodinger equation can be written as 

Let S = r R, then Eq. can be simplified as 



-^S-2fi[XLr-E]S = 0. (59) 



Define a new variable 



A 



^ = (tY [>^Lr-E], (60) 
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such that Eq. ||58| can be transformed as 

S" — xS = 0, 



(61) 



which is the Airy equation. The solution which satisfies the boundary condition S —>■ 
as a; — >^ oo is the Airy function Ai(x). It is easy to show that the eigen-energy formula is 



En Xji [ 



(62) 



where x„ is the n-th zero of the Airy function counting from x = along —x. In refer- 
ence 



IS 



10| , the values Al = 5 and /i = 0.75 are used. In this case, the eigen-energy formula 

^„ = -2.554364772 (63) 



6 Numerical Results 

The matrix elements in the routine dif eq have been calculated in section ^. All but the 
main driver subroutines are left intact. The intial guesses of the wavefunction and its 
derivative in the main driver program are just the square of the sine function and its 
derivative. The codes of the modified driver programs,bohr . c and linear, c, and the 
associating dif eq. c are listed in the appendix of this paper. All the codes are written in 
C. 

For large number of mesh points (1000 < N < 10000), the code generates smooth 
wavefunctions and energies which do not correspond to any eigen modes. An explicit 
Schrodinger solver is known to be conditionally unstable ^. This kind of instabilities 
motivates numerical schemes such as the Numrov method 0, the symplectic scheme 
microgenetic algoritm [^, and others f^. For smaller number of mesh points (40 < N < 
100), it is observed that the relaxed eigenvalues are largely similar to the initial guesses 
and that the smoothness of the relaxed eigenfunctions is very sensitive to the initial 
guesses of the eigenvalue values. Therefore we adopt the criterion that the smoothness 
of the relaxed wavefunction picks out an initial guess as the correct eigenvalue when the 
number of mesh points is small enough. This criterion has certain a priori appeal because 
we expect a real physical wavefunction to be smooth. The reason for using the initial 
guess instead of the relaxed eigenvalue is based on the observation that the former is 
more accurate than the latter when the code is calibrated using exact results. Table 0] 
illustrates the motivation of the smoothness criterion. The combination of a small number 
of mesh points and the smoothness criterion is the simplest way to make the relaxation 
code workable as an explicit Schrodinger solver. 

Figs. [0-0] compare the relaxed and exact wavefunctions of the hydrogen atom for 
the first few quantum numbers. The relaxed wavefunction typically has a vanishingly 
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small amplitude. It is explained by the tendency of the relaxation routine to relax the 
wavefunction to the trivial solution ■0 = 0. In order to compare the relaxed and exact 
wavef unctions, the amplitude of the diminishing relaxed wavefunction is rescaled to match 
the exact wavefunction. Despite the large descrepancies between the relaxed and exact 
wavefunctions of the hydrogen atom in Figs. both the relaxed and exact wavefunc- 
tions share the same basic features. Both relaxed wavefunctions satisfy the smoothness 
criterion. However, there are some difficulties in obtaining a smooth relaxed wavefunc- 
tion for Fig. 0. The relaxed wavefunction in this case also has an unwelcomed kink 
on the right hand side of the plot, which the exact wavefunction does not show. Hence 
the relaxed and exact wavefunctions do not share the same features in this case. The 
discrepancy here is mostly due to degeneracy. The eigen-energy formula of the hydrgon 
atom |]T|] 

Em = ^^jr^ (64) 

reveals the degeneracy over the n and / quantum states. For example, £'21 = E^q such 
that the wavefunctions of these two states are degenerate by virtue of having the same 
eigenvalue. In principle, any linear combination of the degenerate solution vectors also 
satisfies the finite difference equation Eq. |]TU|. It also means that the relaxed wavefunction 



can be any linear combination of the degenerate wavefunctions. This problem will plague 
the numerical solutions of all degenerate states at higher n and I. Imperfect wavefucntions 
is not a liability when the goal is to calculate the eigenvalues in the simplest possible way. 
In this case, the relaxed wavefunction is used merely as a guide to pick out the correct 
eigenvalue and not final product. 

The situation of the numerical solution of the linear potential is slightly better. 
Figs illustrate the relaxed wavefunctions of a linear potential compared to the exact 
wavefunctions. The wavefunctions of a linear potential are not degenerate. The relaxation 
method is expected to be relatively successful in the absence of degeneracy. In this paper, 
we use \l = 5 GeV and n = 0.75 GeV. Table lists the eigenvalues of the first few / 
values in the case of the linear potential. 



7 Conclusion 

In this paper, we show how to circumvent the problem of conditional instability associating 
with an explicit Schrodinger solver when the relaxation method is used. The relaxation 
code with the combination of small number of mesh points and the smoothness criterion 
gives good approximate eigenvalues in both the linear and Coulomb potential cases. The 
role of the relaxed wavefunction is simply to guide the selection of the correct eigenvalue 
and is not used as an answer. More accurate eigenvalues can be obtained by solving the 



momentum space integral equation using the Nystrom method |TT|. In this work, we 



have solved only the non-relativistic Schrodinger equation with the linear and Coulomb 
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potentials. In the case of a relativistic equation, the Hamiltonian contains a kinetic term 
+ — m which is difficult to solve numerically with r-space codes. The p-space 
solution using the Nystrom method on the other hand has the advantages of stability, 
accuracy and robustness over the r-space solution using the relaxation or shooting meth- 
ods. For these reasons, we prefer the p-space codes over the r-space codes for production 
purposes. Nevertheless, the r-space code is useful in checking the p-space results whenever 
the former is available. The numerical r-space calculation is a natural starting point for 
students to learn scientific computation because they would have already seen the exact 
r-space solutions in a traditional quantum mechanics curriculum. Since the numerical 
solution of the hydrogen atom is not readily available in publications, it is hoped that 
the numerical method presented in this paper provides the simplest numerical scheme to 
supplement the discussion of the Schrodinger equation in an advanced undergraduate or 
beginning graduate quantum mechanics course. 
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9 Appendix 



9.1 Non-relativistic Hydrogen Atom 

bohr . c is a driver program for solving the non-relativistic hydrogen atom using the relax- 
ation method. The inputs are n (the principal quantum number), / (the orbital angular 
momentum quantum number), E (the initial guess of the ground state eigen-energy of 
the hydrogen atom which is approximately -13.6) and scale (a scale factor to adjust the 
amplitude of the solution). The output is a file bohr.dat which contains the relaxed 
wavefunction. 

/ * bohr . c */ 

#include <stdio.h> 

#include <math.h> 

#include "difeq.c" 

#include "/recipes/c/solvde.c" 

#include "/recipes/c/bksub . c" 
#include "/recipes/c/pinvs . c" 
#include "/recipes/c/red. c" 
#define NRANSI 

#include "/recipes/c/nrutil .h" 

#include "/recipes/c/nr .h" 

#include "/recipes/c/nrutil . c" 

#define NE 3 

#define M 101 

#define NB 1 

#define NSI NE 

#define NY J NE 

#define NYK M 

#define NCI NE 

#define NCJ (NE-NB+1) 

#define NCK (M+1) 

#define NSJ (2*NE+1) 

#define pi 3.1415926535897932384626433 
int l,inpt=M; 

float h,x[M+l] ,mu,e2,a0; 

int main (void) /* Program bohr.c */ 
{ 
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void solvde(int itmax, float conv, float slowc, float scalv[] , 

int indexv[], int ne, int nb, int m, float **y, float ***c, 
float **s) ; 

int i,itma:x,k,indexv[NE+l] ,n; 

float conv, deriv,sinpx,cospx,ql, slowc, scalv[NE+l] ,eigen, guess; 
float **y,**s,***c, max, scale; 
FILE *fp; 

mu=0.5107208e6; 

e2=7.297353e-3; 

a0=1.0/mu/e2; 

y=matr ix (1 , NY J , 1 , NYK) ; 

s=matrix(l,NSI,l,NSJ) ; 

c=f3tensor(l,NCI,l,NCJ,l,NCK) ; 

itmax=100; 

conv=l . Oe-5 ; 

slowc=1.0; 

h=1.0/(M-l); 

printf("\nn, 1, E, scale = "); 

scanf ("yody„d%fy.f " ,&n,&l,&guess,&scale) ; 

indexv[l]=l; 

indexv [2] =2 ; 

indexv[3] =3; 

for (k=l;k<=(M-l) ;k++) { 
x[k] = (k-l)*h; 
sinpx=sin( (n-1) *pi*x [k] ) ; 
cospx=cos ( (n-1) *pi*x [k] ) ; 
y[l] [k] =sinpx*sinpx; 
y [2] [k] =2 . 0* (n-1) *pi*cospx*sinpx ; 
y[3] [k]=guess/(l+n)/(l+n); 

} 

x[M]=1.0; 
y[l] [M]=0.0; 

y[3] [M]=guess/(l+n)/(l+n) ; 
y[2] [M]=0.0; 
scalv[l]=1.0; 
scalv [2] =1.0; 

scalv [3] =-guess/ (1+n) / (1+n) ; 

solvde (itmajx , conv , slowc , scalv , indexv , NE , NB , M , y , c , s) ; 
eigen=-13 . 598289/ (n+1) / (n+1) ; 

printf ("\nl\tE_in\t\tnumerical E\texact E\t\tnumerical-exact\n") ; 



14 



printf ( "y.d\ty„f \ty„f \ty„f \ty.f \n\n" , 1 , guess/ (1+n) / (1+n) , 

y[3] [1] ,eigen,y[3] [l]-eigeii) ; 
fp=fopen("bohr.dat","w") ; 
max=0 . ; 

for(i=l;i<=M;i++) if (f abs(y [1] [i] )>max) max=f abs (y [1] [i] ) ; 
for(i=l;i<=M;i++) f printf (fp, "y.f\ty.e\n" , (i-l)*h,y[l] [i] /max*scale) ; 
f close(fp) ; 

f ree_f Stensor (c , 1 , NCI , 1 , NCJ , 1 , NCK) ; 
free_inatrix(s,l,NSI,l,NSJ) ; 
f ree.matr ix (y , 1 , NY J , 1 , NYK) ; 
return 0; 

} 

#undef NRANSI 

The dif eq. c code below defines the iS-matrix elements of the finite difference equation 
of the non-relativistic hydrogen atom and is called by the relaxation subroutine solvde . c. 

/* difeq.c for the hydrogen atom driver program */ 

extern int mpt , 1 ; 

extern float h,x[] ,mu,e2,a0; 

void difeq(int k, int kl, int k2, int jsf, int isl, int isf, int indexv[], 
int ne, float **s, float **y) 

{ 

float xk,ylk,y2k,y3k; 
if (k == kl) { 

s[3] [3+indexv[l]]=1.0; 
s[3] [3+indexv[2]]=0.0; 
s [3] [3+indexv [3] ] =0 . ; 
s[3] [jsf]=y[l] [1] ; 
} else if (k > k2) { 

s[l] [3+indexv [1]] =1.0; 

s[l] [3+indexv [2]] =0.0; 

s [1] [3+indexv [3] ] =0 . ; 

s[l] [jsf]=y[l] [mpt] ; 

s[2] [3+indexv [1]] =0.0; 

s[2] [3+indexv [2]] =1.0; 

s [2] [3+indexv [3] ] =0 . ; 

s[2] [jsf]=y[2] [mpt] ; 

} else { 
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s[l] [indexv[l]] = -1.0; 

s[l] [indexv[2]] = -0.5*h; 

s[l] [indexv[3]]=0.0; 

s[l] [3+indexv[l]]=1.0; 

s[l] [3+indexv[2]] = -0.5*h; 

s [1] [3+indexv [3] ] =0 . ; 

xk=(x[k] +x [k-l])/2.0; 

ylk=(y[l] [k]+y[l] [k-l])/2.0; 

y2k= (y [2] [k] +y [2] [k- 1] ) /2 . ; 

y3k=(y[3] [k]+y[3] [k-l])/2.0; 

s [2] [indexv [ 1] ] =h/2 . 0/pow ( 1 . 0-xk ,4.0)* 

(2 . 0*mu*aO*aO* (y3k+ (1 . 0/xk-l . 0) *e2/a0) 
-pow(1.0/xk-1.0,2.0)*l*(l+l)) ; 
s [2] [indexv [2]] = -1 . 0+h/(l .0-xk) ; 
s [2] [indexv [3] ] =h*mu/pow ( 1 . 0-xk ,4.0) *ylk*a0*a0 ; 
s[2] [3+indexv [l]]=s [2] [indexv [1]] ; 
s [2] [3+indexv [2] ] =2 . 0+s [2] [indexv [2] ] ; 
s [2] [3+indexv [3] ] =s [2] [indexv [3] ] ; 
s[3] [indexv [1]] =0.0; 
s[3] [indexv [2]] =0.0; 
s[3] [indexv [3]] = -1.0; 
s [3] [3+indexv [1] ] =0 . ; 
s [3] [3+indexv [2] ] =0 . ; 
s [3] [3+indexv [3] ] =1 . ; 
s[l] [jsf]=y[l] [k]-y[l] [k-l]-h*y2k; 
s[2] [jsf]=y[2] [k]-y[2] [k-1] +2 . 0*h/ (1 . 0-xk) *y2k+ 

h/pow ( 1 . 0-xk , 4 . 0) * (2 . 0*mu*aO*aO* (y3k+ ( 1 . 0/xk- 1 . 0) *e2/a0) 

-pow ( 1 . 0/xk-l . , 2 . 0) *1* (1+1) ) *ylk ; 
s[3] [jsf]=y[3] [k]-y[3] [k-1] ; 
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9.2 NRSE with a Linear Potential 



linear. c is a driver program for solving the non-relativistic Schrodinger equation with 
a hnear potential using the relaxation method. The inputs are n (the principal quantum 
number), / (the orbital angular momentum quantum number), E (the initial guess of the 
eigen-energy) and scale (a scale factor to adjust the amplitude of the wavef unction) The 
output is a file linear.dat which contains the relaxed wavefunction. 

/* linear. c */ 

#include <stdio.h> 
#include <math.h> 
#include "difeq.c" 
#include "/recipes/c/solvde.c" 
#include "/recipes/c/bksub. c" 
#include "/recipes/c/pinvs . c" 
#include "/recipes/c/red. c" 
#define NRANSI 

#include "/recipes/c/nrutil .h" 

#include "/recipes/c/nr .h" 

#include "/recipes/c/nrutil . c" 

#define NE 3 

#define M 101 

#define NB 1 

#define NSI NE 

#define NY J NE 

#define NYK M 

#define NCI NE 

#define NCJ (NE-NB+1) 

#define NCK (M+1) 

#define NSJ (2*NE+1) 

#define pi 3.1415926535897932384626433 
int l,mpt=M; 

float h,x[M+l] ,mu,lambdal; 

int main (void) 
{ 

void solvdeCint itmax, float conv, float slowc, float scalv[] , 

int indexv[] , int ne, int nb, int m, float **y, float ***c, 
float **s) ; 

int i,itma:x,k,indexv[NE+l] ,n; 
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float conv,deriv,sinpx,cospx,ql,slowc,scalv[NE+l] , max, e, scale; 
float **y , **s , ***c ; 
FILE *fp; 

mu=0 . 75 ; 

lambdal=5 . ; 

y=matr ix (1 , NY J , 1 , NYK) ; 

s=matrix(l,NSI,l,NSJ) ; 

c=f3tensor(l,NCI,l,NCJ,l,NCK) ; 

itmax=100; 

conv=l .Oe-6; 

slowc=l . ; 

h=1.0/(M-l); 

printf("\nn, 1, E scale = "); 
scanf ("yod7od%fy.f " ,&n,&l,&e,&scale) ; 
indexv[l]=l; 
indexv [2] =2 ; 
indexv [3] =3 ; 

for (k=l;k<=(M-l) ;k++) { 
x[k] = (k-l)*h; 
sinpx=sin(n*pi*x [k] ) ; 
cospx=cos (n*pi*x [k] ) ; 
y[l] [k] =sinpx*sinpx; 
y[2] [k] =2 . O*n*pi*cospx*sinpx; 
y[3] [k]=e; 

} 

x[M]=1.0; 
y[l] [M]=0.0; 
y[3] [M]=e; 
y[2] [M]=0.0; 
scalv[l]=1.0; 
scalv [2] =1.0; 
scalv[3]=e; 

solvde (itmajx , conv , slowc , scalv , indexv , NE , NB , M , y , c , s) ; 
fp=f open ("linear .dat" , "w") ; 
max=0 . ; 

for(i=l;i<=M;i++) if (fabsCy [1] [i] )>max) max=f abs(y [1] [i] ) ; 

f or (i=l ; i<=M; i++) f printf (f p, '"/.f \t%f \n" , (i-1) *h,y [1] [i] /max*scale) ; 

f close(fp) ; 

printf ("\nl\tE_i\t\tE_f\n") ; 

printf ("y.d\ty.f\ty.f \n\nM,e,y [3] [1]) ; 
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f ree_f atensor (c , 1 , NCI , 1 , NC J , 1 , NCK) ; 
free_matrix(s,l,NSI,l,NSJ) ; 
f ree.matrix (y , 1 , NY J , 1 , NYK) ; 
return 0; 

} 

#uncief NRANSI 

The dif eq. c code below define the S*- matrix elements of the finite difference equation 
of a non-relativistic Schrodinger equation with a linear potential and is called by the 
relaxation subroutine solvde . c. 

/* difeq.c for the linear potential driver program */ 

extern int mpt , 1 ; 

extern float h,x[] ,mu,lainbdal; 

void difeq(int k, int kl, int k2, int jsf, int isl, int isf, int indexvG, 
int ne, float **s, float **y) 

{ 

float xk,ylk,y2k,y3k; 
if (k == kl) { 

s[3] [3+indexv[l]]=1.0; 

s [3] [3+indexv [2] ] =0 . ; 

s [3] [3+indexv [3] ] =0 . ; 

s[3] [jsf]=y[l] [1] ; 
} else if (k > k2) { 



s[l] 


[3+indexv [1]]=1. 


0; 


s[l] 


[3+indexv [2] ] =0 . 


0; 


s[l] 


[3+indexv [3]]=0. 


0; 


s[l] 


[jsf]=y[l] [mpt] ; 




s[2] 


[3+indexv [1]]=0. 


0; 


s[2] 


[3+indexv [2]]=1. 


0; 


s[2] 


[3+indexv [3] ] =0 . 


0; 


s[2j 


[jsf]=y[2] [mpt] ; 




s[l] 


[indexv[l]] = -i 


-.0; 


s[l] 


[indexv[2]] = -0.5*h; 


s[l] 


[indexv[3]]=0.0; 




s[l] 


[3+indexv [1]]=1. 


0; 


s[l] 


[3+indexv [2]] = 


-0.5*h; 


s[l] 


[3+indexv [3] ] =0 . 


0; 
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xk=(x[k] +x [k-l])/2.0; 
ylk=(y[l] [k]+y[l] [k-l])/2.0; 
y2k=(y[2] [k]+y[2] [k-l])/2.0; 
y3k= (y [3] [k] +y [3] [k-1] ) /2 . ; 

s [2] [indexv [1] ] =h/2 . 0/pow(l . 0-xk , 4 . 0) * (2 . 0*mu* 

(y3k+xk/ (xk-1 . 0) *lainbdal) 
-pow(1.0/xk-1.0,2.0)*l*(l+l)) ; 
s [2] [indexv [2]] = -1 .0+h/(l . 0-xk) ; 
s [2] [indexv [3] ] =h*mu/pow ( 1 . 0-xk , 4 . 0) *y Ik ; 
s[2] [3+indexv[l]]=s[2] [indexv [1]] ; 
s [2] [3+indexv [2] ] =2 . 0+s [2] [indexv [2] ] ; 
s [2] [3+indexv [3] ] =s [2] [indexv [3] ] ; 
s[3] [indexv [1]] =0.0; 
s[3] [indexv [2]] =0.0; 
s[3] [indexv [3]] = -1.0; 
s [3] [3+indexv [1] ] =0 . ; 
s [3] [3+indexv [2] ] =0 . ; 
s [3] [3+indexv [3] ] =1 . ; 
s[l] [jsf]=y[l] [k]-y[l] [k-l]-h*y2k; 
s[2] [jsf]=y[2] [k]-y[2] [k-1] +2 . 0*h/(l . 0-xk) *y2k+ 

h/pow ( 1 . 0-xk , 4 . 0) * (2 . 0*mu* (y3k+xk/ (xk- 1.0) *lambdal) 

-powd . 0/xk-l . , 2 . 0) *1* (1+1) ) *ylk ; 
s[3] [jsf]=y[3] [k]-y[3] [k-1] ; 
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X 

Figure 1: The relaxed n = 1, / = wavefunction of the non-relativistic hydrogen atom 
compared against the exact wavefunction. The numerical eigenvalue is —13. 59827 eV, 
versus the exact value of —13.598289 eV. 
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X 

Figure 2: The relaxed n = 2, / = wavefunction of the non-relativistic hydrogen atom 
compared against the exact wavefunction. The numerical eigenvalue is —3. 399750 eV, 
versus the exact value of —3.399572 eV. 
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X 

Figure 3: The relaxed n = 2, / = 1 wavefunction of the non-relativistic hydrogen atom 
compared against the exact wavefunction. The numerical eigenvalue is — 1.510056 eV, 
versus the exact value of — 1.510921 eV. 
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Figure 4: The relaxed / = 0, n = 1 state wavcfunction of a non-rclativistic Schrodinger 
equation with a hnear potential using = 5 GcV and jj, = 0.75 GcV. The numerical 
eigenvalue is 5.9719 GeV, versus the exact value of 5.972379 GeV. The relaxed wavefunc- 
tion has been rescaled to match the exact wavefunction. 
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X 

Figure 5: The relaxed I = 0, n = 2 state wavcfunction of a non-rclativistic Schrodinger 
equation with a hnear potential using = 5 GcV and jj, = 0.75 GeV. The numerical 
eigenvalue is 10.441 GeV, versus the exact value of 10.442114 GeV. The relaxed wave- 
function has been rescaled to match the exact wavef unction. 
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Table 1: Comparisons among the inital guesses, the relaxed and exact eigenvalues of the 
non-rclativistic hydrogen atom and the Schrodinger equation with a linear potential using 
= 5 GeV, n = 0.75 GeV. 









Coulomb 




n 


/ 


Initial 


Numerical 


Exact 


1 





-13.598270 


-13.621142 


-13.598289 


2 





-3.399750 


-3.400535 


-3.399572 


2 


1 


-1.510056 


-1.510060 


-1.510921 








Linear 




n 


/ 


Initial 


Numerical 


Exact 


1 





5.9719 


6.146734 


5.972379 


2 





10.4410 


10.418742 


10.442114 
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Table 2: The numerical eigenvalues of the non-relativistic Schrodinger equation with a 
linear potential using Xl — 5 GeV, /i — 0.75 GeV and n — 1. 



I Numerical 






5.9719 


1 


8.5850 


2 


10.8514 


3 


12.9020 


4 


14.9790 


5 


16.5845 
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